#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include "MSTK.h"


#ifdef __cplusplus
extern "C" {
#endif

int MESH_Surf_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int ring) {
  int iv, i, j, idx, vertex_id, face_id, mkvid, mkfid;
  int nv, nf, nvf, nfv;
  MVertex_ptr mv, mv_new, *fverts, global_mv, local_mv;
  MFace_ptr local_mf, global_mf, mf_new;
  List_ptr global_mvfaces, global_mfvertices,local_mfvertices;
  List_ptr par_vertices, par_faces;
  int *v_g2l, *f_g2l;
  /* basic mesh information */
  /*  RepType rtype = MESH_RepType(mesh);*/
  nv = MESH_Num_Vertices(submesh);
  nf = MESH_Num_Faces(submesh);

  /* store the map between the vertex global ID to ID in submesh */
  v_g2l = (int *)MSTK_malloc((MESH_Num_Vertices(mesh)+1)*sizeof(int));
  /* store the map between the face global ID to ID in submesh */
  f_g2l = (int *)MSTK_malloc((MESH_Num_Faces(mesh)+1)*sizeof(int));

  /* Mark the list of global entities in this submesh*/
  par_vertices=List_New(nv);
  mkvid = MSTK_GetMarker();
  idx = 0;
  while((local_mv = MESH_Next_Vertex(submesh,&idx))) {
    vertex_id = MV_GlobalID(local_mv);
    /* fill the map */
    v_g2l[vertex_id] = MV_ID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    List_Add(par_vertices,global_mv);
    MEnt_Mark(global_mv,mkvid);
  }


  par_faces = List_New(nf);
  mkfid = MSTK_GetMarker();  
  idx = 0;
  while((local_mf = MESH_Next_Face(submesh,&idx))) {
    face_id = MF_GlobalID(local_mf);
    /* fill the map */
    f_g2l[face_id] = MF_ID(local_mf);
    global_mf = MESH_FaceFromID(mesh,face_id);
    List_Add(par_faces,global_mf);
    MEnt_Mark(global_mf,mkfid);
  }

  /*Loop through sub vertices, and see its neighboring faces */
  int index_face = 0;
  for(iv = 0; iv < nv; iv++) {
    /* get global face id and loop through global mesh
     here, can not use List par_vertices, it will be dynamically added */
    local_mv = MESH_Vertex(submesh,iv);
    /* only check vertices on the processor boundary */
    if(MV_PType(local_mv) != PBOUNDARY)
      continue;
    vertex_id = MV_GlobalID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    /* get global neighbor faces of this vertex */
    global_mvfaces = MV_Faces(global_mv);
    nvf = List_Num_Entries(global_mvfaces);
    for(i = 0; i < nvf; i++) {
      global_mf = List_Entry(global_mvfaces,i);
      /* if it is not marked, add this face into submesh */
      if (MEnt_IsMarked(global_mf,mkfid)) {
	local_mf = MESH_FaceFromID(submesh,f_g2l[MF_GlobalID(global_mf)]);
	/* mark PUTILIZED face */
	if(MF_PType(local_mf) != PGHOST) {
	  MF_Set_PType(local_mf,PUTILIZED);
	  local_mfvertices = MF_Vertices(local_mf,1,0);
	  for(j = 0; j < List_Num_Entries(local_mfvertices); j++) {
	    mv = List_Entry(local_mfvertices,j);
	    /* mark PUTILIZED vertex */
	    if(MV_PType(mv) != PBOUNDARY)
	      MV_Set_PType(mv,PUTILIZED);
	  }
	  List_Delete(local_mfvertices);
	}
	continue;
      }
      /* loop through the global vertices of this face */
      global_mfvertices = MF_Vertices(global_mf,1,0);
      nfv = List_Num_Entries(global_mfvertices);
      fverts = (MVertex_ptr *)MSTK_malloc((nfv)*sizeof(MVertex_ptr));
      for(j = 0; j < nfv; j++) {
	global_mv = List_Entry(global_mfvertices,j);
	/* if this vertex is already marked, use its local counterpart */
	if(MEnt_IsMarked(global_mv,mkvid)) {
	  mv = MESH_VertexFromID(submesh,v_g2l[MV_GlobalID(global_mv)]);
	  fverts[j] = mv;
	}
	/* otherwise, add this vertex into submesh */
	if(!MEnt_IsMarked(global_mv,mkvid)) {
	  List_Add(par_vertices,global_mv);
	  MEnt_Mark(global_mv,mkvid);
	  mv_new = MV_New(submesh);
	  PMV_Duplicate(global_mv,mv_new);
	  /* ptype is not copied in PMV_Copy */
	  MV_Set_PType(mv_new,PGHOST);
	  fverts[j] = mv_new;
	  /* fill the map */
	  v_g2l[MV_GlobalID(global_mv)] = MV_ID(mv_new);
	}
      }
      /* add face */
      List_Add(par_faces,global_mf);
      index_face++;
      MEnt_Mark(global_mf,mkfid);
      mf_new = MF_New(submesh);
      MF_Set_Vertices(mf_new,nfv,fverts);
      MF_Set_PType(mf_new,PGHOST);
      PMF_Duplicate(global_mf,mf_new);
      /* fill the map */
      f_g2l[MF_GlobalID(global_mf)] = MF_ID(mf_new);

      if (fverts) MSTK_free(fverts);
      List_Delete(global_mfvertices);
    }
    List_Delete(global_mvfaces);
  }

  List_Unmark(par_vertices,mkvid);
  List_Delete(par_vertices);

  List_Unmark(par_faces,mkfid);
  List_Delete(par_faces);
  MSTK_FreeMarker(mkvid);
  MSTK_FreeMarker(mkfid);
  MSTK_free(v_g2l);  


  return 1;
}
  
#ifdef __cplusplus
}
#endif

